inline	
Riemann_Kernels::Riemann_Kernels(const field<mat> in_riemann_points, int in_dim)
:riemann_points(in_riemann_points), N_points (riemann_points.n_rows), dim(in_dim)
{
  
}

inline
mat
Riemann_Kernels::Stein( std::string tosave_in)
{
  
  cout << "Calculating Stein Kernel " << endl;  
  mat K = zeros(N_points, N_points);
  double sig = 0.5; // Sigma
  double s2 = sig*dim;
  //mat tmp;
  for (uword i = 0; i < N_points ;  ++i)
  {
    for (uword j = i; j < N_points ;  ++j)
    {
      double detX  = pow( det( riemann_points(i,0) ), sig);
      double detY  = pow( det( riemann_points(j,0) ), sig);
      double detXY = pow( det( riemann_points(i,0) + riemann_points(j,0) ), sig);
      K(i,j) = pow(2,dim*sig)*sqrt(detX*detY)/detXY;
      cout << "K(i,j): " << K(i,j) << endl;
      //getchar();
    }
  }
  
  K = symmatu(K);
  K.save(tosave_in, raw_ascii);
  return K;
}


inline
mat
Riemann_Kernels::Stein_faster( std::string tosave_in)
{
  
  cout << "Calculating Stein Kernel Faster?? " << endl;  
  mat K = zeros(N_points, N_points);
  double sig = 0.5; // Sigma
  //double s2 = s*dim;
  double detXpY; //plus
  double detXtY; //times
  double S;
  //mat tmp;
  for (uword i = 0; i < N_points ;  ++i)
  {
    for (uword j = i; j < N_points ;  ++j)
    {
       detXpY = det( (riemann_points(i,0) + riemann_points(j,0))/2 );
       detXtY = det( riemann_points(i,0)*riemann_points(j,0) );

      S = log( detXpY  ) - log(detXtY)/2;
      K(i,j) = exp(-sig*S);
      //cout << "K(i,j): " << K(i,j) << endl;
      //getchar();
    }
  }
  
  K = symmatu(K);
  K.save(tosave_in,raw_ascii);
  return K;
}